Magnetic neutron scattering from spherical nanoparticles with Néel surface anisotropy: atomistic simulations

Based on the Landau–Lifshitz equation, atomistic simulations of the magnetic neutron scattering from inhomogeneously magnetized spherical nanoparticles with a strong surface anisotropy are carried out.


Introduction
Magnetic nanoparticles are the subject of intense worldwide research efforts which are partly motivated by potential applications in areas such as medicine, biology and nanotechnology [see e.g. Lak et al. In the majority of studies, the internal spin structure of the nanoparticles is neglected and assumed to be uniform (called the macro-or superspin model). While this is probably justified in many application-oriented approaches in which an overall understanding is sufficient, it is of interest, at least from the standpoint of fundamental science, to elucidate the effect of a non-uniform spin structure on a certain physical property.
Scattering techniques, in particular employing X-rays and neutrons, have proved to be very powerful in this endeavour, since they provide statistically averaged information on a large number of scattering particles. For instance, using Monte Carlo simulations of a discrete atomistic spin model, Kö hler et al. (2021) have numerically studied the influence of antiphase boundaries in iron oxide nanoparticles on their spin structure. These authors used the Debye scattering equation to relate the internal spin disorder to the broadening of certain X-ray Bragg peaks. Vivas et al. (2020) carried out micromagnetic continuum calculations of the spin structure of defect-free iron nanoparticles and related a vortex-type magnetization configuration to certain signatures in the magnetic neutron scattering cross section and correlation function.
Magnetic small-angle neutron scattering (SANS) is a powerful technique for investigating spin structures on the mesoscopic length scale ($1-100 nm) and inside the volume of magnetic materials Michels, 2021). Recent SANS studies of magnetic nanoparticles, in particular employing spin-polarized neutrons, unanimously demonstrate that their spin textures are highly complex and exhibit a variety of non-uniform, canted or core-shell-type configurations [see e.g. Disch et al. (2012), Krycka et al. (2014),   Honecker et al. (2022) and references therein]. The analysis of magnetic SANS data relies largely on structural form-factor models for the cross section, borrowed from nuclear SANS, which do not properly account for the existing spin inhomogeneity inside a magnetic nanoparticle. Progress in magnetic SANS theory (Honecker & Michels, 2013;Michels et al., 2014;Mettus & Michels, 2015;Erokhin et al., 2015;Metlov & Michels, 2015Michels et al., 2016Michels et al., , 2019Mistonov et al., 2019;Zaporozhets et al., 2022) strongly suggests that, for the analysis of experimental magnetic SANS data, the spatial nanometre-scale variation of the orientation and magnitude of the magnetization vector field must be taken into account, going beyond the macrospin-based models that assume a uniform magnetization.
In this paper, we employ atomistic simulations using the Landau-Lifshitz equation (LLE) to investigate the role of the Né el surface anisotropy in magnetic nanoparticles and its effect on the magnetic SANS cross section and correlation function. We take into account the isotropic exchange interaction, an external magnetic field, a magnetocrystalline anisotropy for the core of the nanoparticles and Né el anisotropy for spins on the surface. The influence of a particle-size distribution function on the magnetic SANS cross section and pair-distance distribution function is also studied. The numerical results reveal marked differences from the superspin model and provide guidance for the experimentalist to identify non-uniform spin structures inside magnetic nanoparticles. We also refer to our accompanying analytical study of the problem , which is restricted to a linear approximation in the magnetization deviation.
The paper is organized as follows. In Section 2 we provide information on the atomistic simulations using the LLE. In Section 3 we display the expressions for the magnetic SANS cross section and for the pair-distance distribution function. The results of the numerical calculations are discussed in Section 4, with Section 4.1 focusing on the effect of the Né el surface anisotropy and Section 4.2 discussing the influence of a log-normal particle-size distribution on the SANS observables. Section 5 summarizes the main findings of this study and provides an outlook on future challenges. Appendix A features results for the SANS observables for different sign combinations of the anisotropy constants.
2. Details of the atomistic SANS modelling using the Landau-Lifshitz equation Fig. 1 shows a schematic depiction of the procedure adopted here to generate and calculate the spin structure, and to obtain the ensuing magnetic SANS cross section and correlation function. This flow-chart-type representation will be discussed in more detail below.
A spherical many-spin nanomagnet is viewed as a crystallite consisting of N atomic magnetic moments l i = a m i , where a denotes the magnitude of the atomic magnetic moment and m i is a unit vector specifying its orientation. We assume the spins 'sit' on a simple cubic lattice, so that a = M s a 3 , where M s is the saturation magnetization of the material and a is the lattice constant. The spherical shape of the nanomagnet is cut from a simple cubic regular grid [ Fig. 1(a)] and its radius R is defined as R ¼ ½ðN À 1Þ=2a, where the integer N is the number of atoms on the side of the cubic grid. The magnetic state of the nanomagnet is investigated with the help of the atomistic approach based on the following Hamiltonian (Dimitrov & Wysin, 1994;Kodama & Berkovitz, 1999;Kachkachi & Garanin, 2001a,b;Iglesias & Labarta, 2001;Kachkachi & Dimian, 2002;Kachkachi & Garanin, 2005;Kazantseva et al., 2008): where H EX is the nearest-neighbour (n.n.) exchange energy, J > 0 is the exchange parameter, H Z denotes the Zeeman energy, B 0 is the homogeneous externally applied magnetic field and H A represents the magnetic anisotropy energy. For the core spins we assume the anisotropy to be of uniaxial symmetry, while for surface spins we adopt the model proposed by Né el (1954). H A;i can then be expressed as where K c > 0 and K s > 0 denote, respectively, the core and surface anisotropy constants, e A is a unit vector along the core anisotropy easy direction, and u ij = (r i À r j )/kr i À r j k is a unit vector connecting the nearest-neighbour spins i and j. The surface spins are defined as those spins which have a coordination number less than six. The magnetodipolar interaction has been ignored in our simulations. This is motivated by the numerical complexity of this energy term, in particular for atomistic simulations (here, for a 10 nm diameter particle the number of spins is N = 11 633), and by the expectation that it is of minor relevance for smaller-sized nanomagnets (Kö hler et al., 2021;Pathak & Hertel, 2021).
The dynamics of each individual magnetic moment m i are described by the Landau-Lifshitz equation (LLE) (Berkov, 2007), where is the gyromagnetic ratio and denotes the damping constant. The deterministic effective magnetic field acting on the spin i is given by The LLE is solved numerically by using the explicit Euler forward-projection method (Baň as, 2005), which consists of two steps. The first step, as seen from equation (6) below, is the simple Euler forward scheme, and the second step, as seen from equation (7), is the projection (or normalization) onto the unit sphere to enforce the constraint km i k = 1. Since we are interested in the static equilibrium, this first-order method is fully appropriate. In equations (6) and (7), k is the time iteration index while i refers to the ith lattice site, h t denotes the time step for the integration procedure. For the termination of the energy minimization, we employ the following criterion: The macroscopic state of the nanomagnet is then described by the following super-or macrospin (representing the net magnetic moment): As an example, we show in Fig. 1(b) the temporal evolution of the Cartesian magnetization components of m and in Fig. 1(c) the numerically computed equilibrium spin configuration for a spherical nanomagnet at zero applied field in a plane across its centre. It is seen that the spins at the centre of the nanoparticle are directed along m while the surface spins exhibit significant misalignment, which is due to the presence of the Né el surface anisotropy. Note that m i are unit vectors, whereas generally kmk 6 ¼ 1.
In our simulations we use the following parameters: atomic magnetic moment a = 1.577 Â 10 À23 A m 2 (corresponding to 1.7 B with B the Bohr magneton), lattice constant a = 0.3554 nm, M s = a /a 3 = 351 kA m À1 , exchange constant J = 8.7 Â 10 À22 J atom À1 , core anisotropy constant K c = 3 Â 10 À24 J atom À1 , damping constant = 3 Â 10 11 s À1 T À1 , gyromagnetic constant = 1.76 Â 10 11 s À1 T À1 and an integration time step of h t = 5 fs. The surface anisotropy constant K s was   used as an adjustable parameter. Experimental K s values for nanoparticles and thin films can be found in the work of Gradmann (1986), Batlle et al. (2022) andO'Handley (2000). A value of K s = 5.22 Â 10 À21 J atom À1 has been estimated by Kachkachi & Dimian (2002) for a 4 nm-sized face-centred cubic cobalt particle.
For the calculation of the magnetic SANS cross section dAE M /d [ Fig. 1( f )], it is necessary to compute the discrete Fourier transform of all the m i belonging to the spherical nanomagnet [ Fig. 1(e)]. In Section 3, the expressions for dAE M / d are formulated for a continuous magnetization distribution M(r) and of its Fourier transform e M MðqÞ. These functions are defined as follows: Using where r i is the location point of the ith spin and q represents the wavevector (scattering vector, defined in Fig. 2). Equation (12) establishes the relation between the outcome of the simulations, m i , and the magnetic SANS cross section, dAE M / d. In the standard SANS geometry, the q space of interest is defined by q = q[0, sin , cos ], which corresponds to the two-dimensional detector plane (q x = 0, see Fig. 2 At each value of the external field, atomistic simulations of the spin structure and of the ensuing magnetic SANS cross section were carried out for 256 random orientations of the core anisotropy axes e A of the particle with respect to the field B 0 . More specifically, once the lattice orientation has been randomly selected, the easy-axis orientation of the particle's core and the distribution of the Né el anisotropy are fixed. The whole system (core plus surface anisotropy) is then randomly rotated relative to B 0 . For the generation of the random angles, we used the low-discrepancy Sobol sequence (sob, https://www.mathworks.com/help/stats/sobolset.html). Therefore, except Fig. 3, all the data shown in this paper correspond to an ensemble of randomly oriented particles. The simulations were carried out by starting from a large positive (saturating) field of about 10 T, and then the field was reduced in steps of, typically, 30 mT.

Magnetic SANS cross section and pair-distance distribution function
The quantity of interest in experimental SANS studies is the elastic magnetic differential scattering cross section dAE M /d, which is usually recorded on a two-dimensional positionsensitive detector. For the most commonly used scattering geometry in magnetic SANS experiments, where the applied magnetic field B 0 k e z is perpendicular to the wavevector k 0 k e x of the incident neutrons (see Fig. 2 where V is the scattering volume, b H = 2.91 Â 10 8 A À1 m À1 is the magnetic scattering length in the small-angle regime (the atomic magnetic form factor is approximated by 1 since we are dealing with forward scattering), e M MðqÞ = ½ e M M x ðqÞ; e M M y ðqÞ; e M M z ðqÞ represents the Fourier transform of the magnetization vector field M(r) = [M x (r), M y (r), M z (r)], denotes the angle between q and B 0 , and the asterisk * stands for the complexconjugated quantity. Note that in the perpendicular scattering geometry the Fourier components are evaluated in the plane q x = 0 (see Fig. 2).
The numerically computed magnetic SANS cross sections that are displayed in this paper correspond to the following average: A sketch of the neutron scattering geometry. The applied magnetic field B 0 k e z is perpendicular to the wavevector k 0 k e x of the incident neutron beam (B 0 ? k 0 ). The momentum transfer or scattering vector q is defined as the difference between k 0 and k 1 , i.e. q = k 0 À k 1 . SANS is usually implemented as elastic scattering (k 0 = k 1 = 2/) and the component of q along the incident neutron beam, here q x , is much smaller than the other two components, so that q ffi ½0; q y ; q z ¼ q½0; sin ; cos . This demonstrates that SANS predominantly probes correlations in the plane perpendicular to the incident beam. For elastic scattering, the magnitude of q is given by q ¼ ð4=Þ sinð Þ, where denotes the mean wavelength of the neutrons and 2 is the scattering angle. The angle = /(q, B 0 ) is used to describe the angular anisotropy of the recorded scattering pattern on the two-dimensional position-sensitive detector.
where dAE M, k /d represents (for fixed K s and B 0 ) the magnetic SANS cross section for a particular core easy-axis orientation e A (with reference to the index 'k') and K denotes the number of random configurations. Equation (14) implies the absence of interparticle interactions. For a uniformly magnetized spherical particle with its saturation direction parallel to e z , i.e. M x = M y = 0 and M z = M s , equation (13) where V p = 4R 3 / 3 is the particle's volume, ðÁÞ 2 mag = b 2 H ðÁMÞ 2 = b 2 H M 2 s is the magnetic scattering-length density contrast and j 1 (qR) is the first-order spherical Bessel function. The well known analytical result for the homogeneous sphere case [equation (15)] and its correlation function [see equation (18) below] serve as a reference for comparison with the nonuniform case.
It is often convenient to average the two-dimensional SANS cross section (dAE M /d)(q) = (dAE M /d)(q y , q z ) = (dAE M /d)(q, ) along certain directions in q space, e.g. parallel ( = 0) or perpendicular ( = /2) to the applied magnetic field, or even over the full angular range. In the following, we consider the 2 azimuthally averaged magnetic SANS cross section, which is used to compute the pair-distance distribution function p(r) according to where j 0 ðqrÞ ¼ sinðqrÞ=ðqrÞ is the spherical Bessel function of zero order. p(r) corresponds to the distribution of real-space distances between volume elements inside the particle weighted by the excess scattering-length density distribution; see the reviews by Glatter (1982) and Svergun & Koch (2003) for detailed discussions of the properties of p(r) and for information on how to compute it by indirect Fourier transformation (Bender et al., 2017). For our discrete simulation data, the integrals in equations (16) and (17) were approximated by the trapezoidal rule. Apart from constant prefactors, p(r) of the azimuthally averaged single-particle cross section [equation (15)], corresponding to a uniform sphere magnetization, is given by (for r 2R) We also display results for the correlation function c(r), which is related to p(r) by As we will demonstrate in the following, when the particles' spin structure is inhomogeneous, dAE M /d and the corresponding p(r) and c(r) differ significantly from the homogeneous case [equations (15) and (18)], which serves as a reference. Because of the r 2 factor, features in p(r) at medium and large distances are more pronounced than those in c(r).

Results and discussion
4.1. Effect of the Néel surface anisotropy The computed normalized magnetization m z [compare equation (9)] of an ensemble of randomly oriented spherical nanomagnets for different values of the surface anisotropy constant K s (in units of 10 À23 J atom À1 , see inset). The particle diameter is D = 10 nm and the remanence values are indicated in the inset. Selected 3D equilibrium spin structures arising from the Né el surface anisotropy [compare also with Figs. 2 and 3 in the accompanying analytical study ]. (a) K s = 5.22 Â 10 À23 J atom À1 and (b) K s = 52.2 Â 10 À23 J atom À1 . Further parameters are core-anisotropy axis e A = [0, 0, 1], core-anisotropy constant K c = 3 Â 10 À24 J atom À1 and external magnetic field B 0 = [0, 0, 150 mT]. The particle diameter is D = 5 nm. The colour code depicts the spin misalignment relative to the average magnetization vector, namely m j = km j À m=kmkk. At the surface of the nanomagnet the spin deviations are larger than those in the core.
surface anisotropy constant K s , and Fig. 4 shows computed hysteresis curves for an ensemble of randomly oriented 10 nmsized nanomagnets. As expected, increasing K s results, for a given particle size, in a progressive surface spin disorder which propagates into the bulk of the nanomagnet. The effect of an enhanced K s also becomes visible in the magnetization curves via an increased coercivity H c and remanence m r . For K s = 0 and dominant exchange, we recover the well known results from the Stoner-Wohlfarth model (Usov & Peschany, 1997), i.e. we find a reduced remanence of m r = 0.5 and a coercivity of where N c denotes the number of atoms belonging to the particle's core. Note that for the case of a strong surface anisotropy [ Fig. 3(b)], the mean magnetization at remanence deviates strongly from the core anisotropy axis (parallel to e z ), which is in contrast to the case of weak anisotropy [ Fig. 3(a)].
This observation is in agreement with the analytical calculations by Garanin & Kachkachi (2003) who predicted the emergence of an effective anisotropy of cubic symmetry for dominant K s . Therefore, with increasing K s we initially observe in Fig. 4 an increase in the remanence. However, for the largest K s , the reduced remanence again decreases slightly from 0.72 to 0.70. We believe that this observation is due to the disordering effect of the surface anisotropy beyond a certain critical K s .  (13)]. Fig. 6 shows the corresponding plots at a (nearly) saturating field of B 0 = 10 T. We emphasize that the depicted scalar functions represent projections of the corresponding three-dimensional quantities onto the q y q z detector plane at q x = 0 (see Fig. 2 (21) and associated discussion in the main text]. The CT (and hence the corresponding ) can take on negative values, but in this figure we show (due to the chosen logarithmic colour scale) the absolute value of the CT. The data correspond to an ensemble of randomly oriented 10 nm-sized nanomagnets. The K s values for each row are (first row) K s = 0, (second row) K s = 5.22 Â 10 À23 J atom À1 , (third row) K s = 26.1 Â 10 À23 J atom À1 and (fourth row) K s = 52.2 Â 10 À23 J atom À1 . anisotropy constant K s increases from the top to the bottom row in Figs. 5 and 6. It is seen that, generally, all the Fourier components contribute to dAE M /d.
Near saturation (Fig. 6), dAE M /d is dominated for all values of K s by the isotropic (-independent) j e M M z j 2 Fourier component and exhibits the characteristic sin 2 anisotropy with two maxima along the vertical direction [compare equation (13)]. Increasing K s enhances the contributions of both transverse Fourier components j e M M x j 2 and j e M M y j 2 and of the CT. Moreover, the latter contributions develop a pronounced angular anisotropy with increasing K s .
At remanence (Fig. 5), dAE M /d and all the Fourier components are isotropic for small values of K s and become progressively more anisotropic with increasing K s . For instance, j e M M z j 2 is initially isotropic and develops a pronounced angular anisotropy that is elongated along the q y direction for larger K s . The CT also develops an anisotropy with increasing K s , with maxima roughly along the detector diagonals. An anisotropic magnetic SANS cross section at zero applied magnetic field of an ensemble of randomly oriented nanoparticles has also been found in the micromagnetic continuum simulations of Vivas et al. (2020). These authors did not consider the Né el surface anisotropy but included the magnetodipolar interaction.
To quantify the fraction of the individual Fourier components in equation (13) relative to the total magnetic SANS cross section dAE M /d, we compute the following dimensionless quantity: where (q, ) is, respectively, given by Kj e M M x j 2 , Kj e M M y j 2 cos 2 , Kj e M M z j 2 sin 2 and K CT sin cos , with K ¼ 8 3 b 2 H V À1 . q max is taken as 10 nm À1 . The corresponding numbers are given as % values in Figs. 5 and 6, and we note that the contribution related to K CT sin cos can be positive as well as negative, in contrast to the other three contributions which are strictly positive. Using the inequality j e M M y cos À e M M z sin j 2 ! 0, it can easily be shown that the contribution K CT sin cos is, however, always smaller than the sum of the other terms (as it must be). We emphasize that the colour-coded plots in Figs. 5 and 6 show the respective Fourier components without the trigonometric functions in equation (13), whereas the quantities do contain the trigonometric terms. For K s = 0 and zero field, the contributions of j e M M x j 2 , j e M M y j 2 and j e M M z j 2 to dAE M / d are approximately equal (while CT = 0). This can be understood by noting the isotropy of these functions and by research papers taking into account the trigonometric terms cos 2 (for j e M M y j 2 ) and sin 2 (for j e M M z j 2 ), which yield a factor of 1/2 on azimuthal averaging [ integration, compare equation (21)]. At remanence, the results for the smallest nonzero K s are nearly identical to the data for K s = 0 (compare Fig. 4).
The effect of increasing K s on the 2 azimuthally averaged I(q) = dAE M /d and on p(r) and c(r) is shown in Fig. 7 for the remanent state and in Fig. 8 for B 0 = 10 T. With increasing spin disorder (induced by an increasing K s ) we observe in Fig. 7(a) that (i) the characteristic form-factor oscillations of I(q) are progressively damped and (ii) the extrema in I(q) shift to larger q values due to the reduced coherent magnetic size of the particle. Generally, in experimental situations, the smearing of form-factor oscillations is related to the effect of a particle-size distribution function and/or experimental resolution. Therefore, when data such as those in Fig. 7(a) are fitted to a set of single-domain particles with a distribution of sizes, rather than to a set of non-uniformly magnetized particles that all have the same size, an erroneous value for the particle size may result. At (quasi)saturation [ Fig. 8(c)] and for small K s at remanence [ Fig. 7(c)], we recover the analytically known expressions for I(q), p(r) and c(r) for uniformly magnetized spherical particles [equations (15) and (18) The effect of the surface anisotropy constant K s (in units of 10 À23 J atom À1 , see inset) on (a) the azimuthally averaged magnetic SANS cross section I(q) = (dAE M /d)(q) (log-log scale), (b) the value of the magnetic SANS cross section at the origin, I(q = 0) versus K s , (c) the pair-distance distribution function p(r) and (d) the correlation function c(r). The data correspond to the remanent state (B 0 = 0 T) and the nanomagnets' diameter is 10 nm. The green dashed line in panel (c) displays the analytical pair-distance distribution function for the case of a uniformly magnetized spherical particle [proportional to equation (18)], where the magnitude is normalized to the maximum value from the numerical simulation in the case K s = 0.

Figure 8
The same as Fig. 7, but for B 0 = 10 T. magnetization per atom (Marshall & Lowde, 1968). We see that, as expected, the increase in K s has a large effect on (0), whereas the reduction is relatively small at 10 T. We also refer to Fig. 10 in Appendix A, where the results for the SANS observables are shown for the other possible sign combinations of the core and surface anisotropy constants.

Effect of a particle-size distribution
In SANS experiments on nanoparticles one always has to deal with a distribution of particle sizes and shapes. The size of a particle has an important effect on its spin structure, e.g. smaller particles generally tend to be nearly uniformly magnetized (due to the dominant role of the exchange interaction), whereas larger particles may exhibit highly inhomogeneous spin structures (due to the magnetodipolar interaction) (Vivas et al., 2020). It is therefore also of interest to study the influence of a distribution of particle sizes on the magnetic SANS observables [dAE M /d, p(r), c(r)]. This has been done using a log-normal probability distribution function, which is defined as (Krill & Birringer, 1998) where denotes the expectation value and 2 is the variance, such that where the corresponding median Ã is determined by the relation For given values of and , the average magnetic SANS cross section h . . . i is computed as P ' denotes the probability related to the particle-size class D ' = 2R ' (diameter), which is computed as The effect of a log-normal particle-size distribution function on the SANS observables (K s = 52.2 Â 10 À23 J atom À1 ). Shown are the two-dimensional hdAE M /di, the corresponding azimuthally averaged hI(q)i, the pair-distance distribution functions hp(r)i and the particle-size distributions w(D) for values of (a) 3 nm, (b) 2 nm and (c) 1 nm. These values correspond to [ln (1 + 2 / 2 )] 1/2 values of, respectively, 0.29, 0.20 and 0.10. The nanoparticles' mean diameter (expectation value) was chosen as = 10 nm in each case. The data correspond to the remanent (B 0 = 0 T) and saturated (B 0 = 10 T) magnetization states. The discrete particle-size classes are defined by the particle diameters D = 3-16 nm with an equidistant step size of ÁD = 1 nm. The black dashed hp(r)i curves are the analytically known solutions for uniformly magnetized spheres of size = 10 nm in the fully saturated state. Fig. 9 summarizes the results obtained for the magnetic SANS cross section and correlation function. As expected, one observes a smearing of the SANS signal with increasing standard deviation of the distribution, which becomes particularly visible in the azimuthally averaged hI(q)i curves via the suppression of the form-factor oscillations. The angular anisotropy of the SANS cross section in the remanent state, which can be seen as a characteristic signature of the Né el surface anisotropy (compare also the lowest row in Fig. 5), becomes less pronounced for large . An increasing applied field suppresses the internal spin disorder and in this way increases the coherent magnetic sizes of the nanoparticles, so that the maximum of hp(r)i shifts to larger distances. At the same time, an increasing field also suppresses fluctuations in the local magnetizations relative to the mean directions, which then results in a reduction in the magnitude of hp(r)i. Up to this point, we have exclusively considered the case where both anisotropy constants are positive, i.e. K c > 0 and K s > 0 (for fixed magnitudes). In Appendix A we also show the results for the SANS observables for the other possible sign combinations of the anisotropy constants. There, we can see that the angular anisotropy and the q dependence of hdAE M / di may be taken as an indication for distinguishing between the cases of positive and negative K s .
In contrast to our accompanying analytical study (Adams et al., 2022), which is based on the linearization of Brown's static equations of micromagnetics, the present numerical work takes into account the full nonlinearity of the underlying equations for the spin dynamics via the Landau-Lifshitz equation. Therefore, the analytical approach is only valid for weak surface anisotropy, while the numerical approach considers surface anisotropy of arbitrary strengths. In both studies, the dipolar interaction is neglected, which is related to its mathematical complexity and the enormous numerical effort to include it in atomistic simulations of large particles.

Conclusions and outlook
We have studied the spin structure and magnetic neutron scattering signal of an ensemble of randomly oriented spherical nanomagnets using the Landau-Lifshitz equation, with particular focus on the Né el surface anisotropy. Taking   account the isotropic exchange interaction, an external magnetic field, a uniaxial magnetic core anisotropy and the Né el surface anisotropy, we compute the magnetic small-angle neutron scattering cross section and the pair-distance distribution function from the obtained equilibrium spin structures. The numerical results are compared with the well known analytical expressions for uniformly magnetized particles. With increasing internal spin disorder (increasing surface anisotropy K s ), the pair-distance distribution function (at remanence) exhibits a systematic shift of its maximum to smaller r values and the total magnetic SANS cross section develops a characteristic anisotropic scattering pattern. The strength of the simulation methodology is that the field evolution of the individual Fourier components and their contribution to the magnetic SANS signal can be monitored. Atomistic and micromagnetic continuum simulations have contributed and will continue to contribute to the fundamental understanding of magnetic SANS.
In our future work, we will focus on the inclusion of both the intraparticle and interparticle dipole-dipole energy and the Dzyaloshinskii-Moriya interaction, which will give rise to more complicated spin textures (e.g. vortex-type structures), in particular for larger particle sizes. Moreover, it is of interest to compare the Né el anisotropy with other phenomenological expressions for the surface anisotropy, such as energy densities of the type AE 1 2 K s ðm Á nÞ 2 , where n is the normal unit vector to the surface (instead of u ij ), or with the case of a truly random surface anisotropy, where u ij are random vectors. In this regard, the present first atomistic simulations may be considered as the starting point towards a more complete description of magnetic SANS.
The supporting information to this paper features a video that displays the SANS observables during the magnetizationreversal process for the case of a strong surface anisotropy.

APPENDIX A Summary of results for different signs of the anisotropy constants
In the main text, we have exclusively considered the case that both anisotropy constants are positive, i.e. K c > 0 and K s > 0. Here, for completeness, we also display in Fig. 10 the results for the SANS observables for the cases of K c > 0 and K s < 0, K c < 0 and K s > 0, and K c < 0 and K s < 0. By comparison with equation (3) it is seen that changing the sign of the core anisotropy constant K c from positive to negative changes the orientation from easy axis to easy plane. Likewise, changing the sign of K s from positive to negative favours the m i k u ij orientation (which we term the tangential orientation) over the m i ? u ij (radial) orientation. The magnitudes of the anisotropy constants are |K c | = 3 Â 10 À24 J atom À1 and |K s | = 52.2 Â 10 À23 J atom À1 . A log-normal particle-size distribution function with a mean diameter of = 10 nm and a width of = 1 nm has been assumed.
It is seen in Fig. 10 that the sign change of K s has the largest effect on the SANS observables. The surface spin structure [ Fig. 10(a)] changes from a more radial spin orientation for K s > 0 to a more tangential orientation for K s < 0. The corresponding angular anisotropy of the two-dimensional hdAE M /di [ Fig. 10(b)] changes from a horizontally elongated (K s > 0) to a vertically elongated pattern (K s < 0). The sign change of K s also manifests in the hp(r)i and hc(r)i functions, but the most prominent effect is seen in hI(q)i [ Fig. 10(c)], where at q ffi 0.8-1.4 nm À1 the functional dependency of hI(q)i is distinctly different, more shoulder-like for K s > 0 to more peak-like for K s < 0 [see inset in Fig. 10(c)]. This feature could be taken as an indication for distinguishing between the two cases of K s > 0 and K s < 0.

Funding information
Funding for this research was provided by National Research Fund of Luxembourg (grant No. 15639149 to MPA and AM).